/*
 * OptionPricer.cpp
 *
 *  Created on: Nov 12, 2009
 *      Author: karthik
 */
#include <fstream>
#include "OptionPricer.h"
#include "halton.h"
#include "mtrand.h"
#include "Date.h"
#include "normal.h"
#include "Utils.h"

double FinalValue(double spot, double volatility, double rate, Date current_date, Date expiration, double randx){
	//	std::cout<<"Inside BS:"<<endl;
	double final_value;
	double tau = time_diff_yrs(expiration,current_date);
	final_value = spot*exp( ((rate-(pow(volatility,2)/2))*(tau)) + volatility*sqrt(tau)*randx );
	//	std::cout<<"final_value:"<<final_value<<endl;
	return final_value;
}


OptionPricer::OptionPricer(){

}
/*
OptionPricer::OptionPricer(Option opt,vector<double> starting_prices,double rate,vector<double> vols,Array2D<double> corr_matrix){

}*/

void OptionPricer::initialize(Option opt,vector<double> starting_prices,double rate,vector<double> vols,Array2D<double> corr_matrix){
	this->opt=opt;
	this->starting_prices=starting_prices;
	this->rate=rate;
	this->vols=vols;
	this->corr_matrix=corr_matrix;
}

double OptionPricer::price(string randGenMethod, int trials, Date present_date){

	double ret=0.0;
	int n=starting_prices.size();	// no of stocks in the basket.


	string filename ="Running_averages_"+randGenMethod+".csv";
	ofstream run_avgs(filename.c_str());

	//	Cholesky Decomposition.



//	cout<<"Cholesky :::::::::::"<<endl;


	Eigenvalue<double> E(corr_matrix);

	Array1D<double> d(n,0.0);
	Array2D<double> U(n,n,0.0);
	Array2D<double> D(n,n,0.0);
	Array2D<double> Dk(n,n,0.0);
	Array2D<double> Ut(n,n,0.0);
	Array2D<double> Ck(n,n,0.0);
	Array2D<double> corr_matrix_new(n,n,0.0);

	E.getRealEigenvalues(d);
	E.getV(U);
	E.getD(D);




	Array1D<double> d_cum(n,0.0);
	Array1D<double> d_new(n,0.0);	//	has eigen values for the first k elements whose sum > 75% of n.


	double sum=0.0;
	for(int i=d.dim()-1;i>=0;i--){
		//		cout << " : " << d[i] ;
		sum+=d[i] ;
		d_cum[i]=sum;
		d_new[i]=d[i];
//		cout << " : " << d_new[i] ;
		Dk[i][i]=D[i][i];
		if (sum>=0.75*n){

			break;
		}
	}
	cout << endl;

	for(int i=0;i<n;i++){
		for(int j=i;j<n;j++){
			Ut[i][j]=U[j][i];
			Ut[j][i]=U[i][j];
		}
	}


	Ck = matmult(matmult(U,Dk),Ut);
	corr_matrix_new = matmult(matmult(U,D),Ut);

	cout<<"\nCorr matrix :::::::::::"<<endl;
	printArray2D(corr_matrix);

	//	printArray1D(d);
	//	printArray2D(U);

	cout<<"\nD (complete) :::::::::::"<<endl;
	printArray2D(D);

//	printArray2D(Ut);
//	cout<<"C(new) :::::::::::"<<endl;
//	printArray2D(corr_matrix_new);
//	cout<<"Ck :::::::::::"<<endl;
//	printArray2D(Ck);

	Cholesky<double> C(corr_matrix_new);

	Array2D<double> L(n,n,0.0);
	L= C.getL();


	PayOff p1;
	p1.strike = opt.strike;

	Halton H(n);
	H.getNext();	// to avoid first vector with all zeroes, which gives nan when taken inverse !!!

	vector<double> running_avg;


	/*
	 * If running for uncorrelated increments, uncomment this part.
	 *
	L = Array2D<double>(n,n,0.0);
	for(int i=0;i<n;i++){
			L[i][i]=1.0;
	}
	*/

	cout << "\nChol :" << L << endl;

	double sum_of_payoffs=0;
	for (int j=0;j<trials;j++){
//				std::cout<<"Trial:"<<j<<endl;

		double corrFinalPrices[n] ;
		vector<double> myRands;

		if (randGenMethod!="Halton"){
			MTRand drand; // double in [0, 1) generator, already init
			vector<double> mtRands;
			for(int i=0;i<n;i++){
				mtRands.push_back(normsinv(drand()));
				//			std::printf("%10.8f ", mtRands[i]);
			}
			//			cout<<"Using mersenne-twister :";
			//			printVofDouble(mtRands);

			myRands = mtRands;
		}else{
			//		Array1D<double> hRands(n,0.0);
			vector<double> hRands=H.getNextGaussian();
			//vector<double> hRands=normsinv_vector(H.getNext());
			//			cout<<"Using halton :";
//						printVofDouble(hRands);

			myRands = hRands;
		}


		vector<double> corrRands ;
		for (int k=0;k<n;k++){
			double rand=0.0;
			for (int l=0;l<n;l++){
				rand+=L[k][l]*myRands[l];
			}
			corrRands.push_back(rand);
		}

		for (int i = 0; i < n; ++i) {
			corrFinalPrices[i]=FinalValue(starting_prices[i], vols[i], rate, present_date, opt.expiration, corrRands[i]);
		}
		//double payoff = p1.getPayOff(corrFinalPrices,n);
		//cout<<"Payoff : "<<payoff<<endl;
		sum_of_payoffs += p1.getPayOff(corrFinalPrices,n);
		running_avg.push_back(sum_of_payoffs/(j+1));

		if (j%100==0){
			run_avgs<<j<<","<<sum_of_payoffs/(j+1)<<"\n";
		}

	}

	run_avgs.close();


	double tau = time_diff_yrs(opt.expiration,present_date);
	cout<<"Expiration time (yrs) : "<<tau<<endl;
	ret = (sum_of_payoffs/trials)*exp(-1*rate*tau);
//	cout<<"Multifactor Option Price : "<<ret<<endl;

	return ret;
}

OptionPricer::~OptionPricer(){

}
